require(plm)
require(runjags)
require(coda)
require(wfe)

#### data preparation ####
LDP.data.Study1 <- read.csv("LDP_data_Study1.csv")

# exclude periods with no speech
# PS = plenary session, BC = Budget Committee, OC = other committees
LDP.personnel.period.date.Study1 <- as.Date(c("1979-11-09", "1980-07-17", "1981-11-30", 
                                              "1982-11-27", "1983-12-27", "1984-11-01", 
                                              "1985-12-28", "1986-07-22", "1987-11-06", 
                                              "1988-12-27", "1989-06-03", "1989-08-10", 
                                              "1990-02-28", "1990-12-29", "1991-11-05", 
                                              "1992-12-12", "1993-08-09", "1993-09-30", 
                                              "1994-07-05", "1995-01-19", "1995-10-01", 
                                              "1996-01-11", "1996-11-07", "1997-09-11", 
                                              "1998-07-30", "1999-01-14", "1999-10-05", 
                                              "2000-04-05", "2000-07-04", "2000-12-05", 
                                              "2001-04-26", "2002-09-30", "2003-09-22", 
                                              "2003-11-19", "2004-09-27", "2005-09-21", 
                                              "2005-10-31", "2006-09-26"))

LDP.PS.data.Study1 <- LDP.BC.data.Study1 <- LDP.OC.data.Study1 <- LDP.data.Study1
for (i in 1:(length(LDP.personnel.period.date.Study1) - 1)) {
  temp.data <- subset(LDP.data.Study1, time == i)
  if (sum(temp.data$PS.char > 0) == 0) {
    LDP.PS.data.Study1 <- subset(LDP.PS.data.Study1, time != i)
  }
  if (sum(temp.data$BC.char > 0) == 0) {
    LDP.BC.data.Study1 <- subset(LDP.BC.data.Study1, time != i)
  }
  if (sum(temp.data$OC.char > 0) == 0) {
    LDP.OC.data.Study1 <- subset(LDP.OC.data.Study1, time != i)
  }
}

# split data by electoral systems
LDP.PS.SNTV.data.Study1 <- subset(LDP.PS.data.Study1, MMM == 0)
LDP.PS.MMM.data.Study1 <- subset(LDP.PS.data.Study1, MMM == 1)
LDP.BC.SNTV.data.Study1 <- subset(LDP.BC.data.Study1, MMM == 0)
LDP.BC.MMM.data.Study1 <- subset(LDP.BC.data.Study1, MMM == 1)
LDP.OC.SNTV.data.Study1 <- subset(LDP.OC.data.Study1, MMM == 0)
LDP.OC.MMM.data.Study1 <- subset(LDP.OC.data.Study1, MMM == 1)

# convert data frames to the pdata.frame class of the plm package
LDP.PS.SNTV.pdata.Study1 <- pdata.frame(LDP.PS.SNTV.data.Study1, c("name_jp", "time"))
LDP.BC.SNTV.pdata.Study1 <- pdata.frame(LDP.BC.SNTV.data.Study1, c("name_jp", "time"))
LDP.OC.SNTV.pdata.Study1 <- pdata.frame(LDP.OC.SNTV.data.Study1, c("name_jp", "time"))
LDP.PS.MMM.pdata.Study1 <- pdata.frame(LDP.PS.MMM.data.Study1, c("name_jp", "time"))
LDP.BC.MMM.pdata.Study1 <- pdata.frame(LDP.BC.MMM.data.Study1, c("name_jp", "time"))
LDP.OC.MMM.pdata.Study1 <- pdata.frame(LDP.OC.MMM.data.Study1, c("name_jp", "time"))

#### main analyses (Table 2) ####
## SNTV period
Study1.LDP.PS.SNTV <- plm(PS.char ~ leader + chair + minister + 
                            age + log.term + vote.share + 
                            E.37 + E.38 + E.39 + E.40, 
                          data = LDP.PS.SNTV.pdata.Study1, 
                          model = "within")
round(summary(Study1.LDP.PS.SNTV, 
              vcov = vcovHC(Study1.LDP.PS.SNTV, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)
nrow(LDP.PS.SNTV.pdata.Study1)  # number of observations
length(unique(LDP.PS.SNTV.pdata.Study1$name_jp))  # number of MPs
length(unique(LDP.PS.SNTV.pdata.Study1$time))  # number of periods

round(exp(Study1.LDP.PS.SNTV$coefficients[1]), 3)  # substantive interpretation

## MMM period
Study1.LDP.PS.MMM <- plm(PS.char ~ leader + chair + 
                           minister + jr.minister + age + 
                           log.term + vote.share * SMD + zombie + 
                           E.42 + E.43 + E.44, 
                         data = LDP.PS.MMM.pdata.Study1, 
                         model = "within")
round(summary(Study1.LDP.PS.MMM, 
              vcov = vcovHC(Study1.LDP.PS.MMM, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)
nrow(LDP.PS.MMM.pdata.Study1)  # number of observations
length(unique(LDP.PS.MMM.pdata.Study1$name_jp))  # number of MPs
length(unique(LDP.PS.MMM.pdata.Study1$time))  # number of periods

round(exp(Study1.LDP.PS.MMM$coefficients[1]), 3)  # substantive interpretation

#### analysis of committees (Online Appendix C) ####
## SNTV, Budget Committee
Study1.LDP.BC.SNTV <- plm(BC.char ~ leader + chair + minister + 
                            age + log.term + vote.share + 
                            E.37 + E.38 + E.39 + E.40, 
                          data = LDP.BC.SNTV.pdata.Study1, 
                          model = "within")
round(summary(Study1.LDP.BC.SNTV, 
              vcov = vcovHC(Study1.LDP.BC.SNTV, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)
nrow(LDP.BC.SNTV.pdata.Study1)  # number of observations
length(unique(LDP.BC.SNTV.pdata.Study1$name_jp))  # number of MPs
length(unique(LDP.BC.SNTV.pdata.Study1$time))  # number of periods

## SNTV, other committees
Study1.LDP.OC.SNTV <- plm(OC.char ~ leader + chair + minister + 
                            age + log.term + vote.share + 
                            E.37 + E.38 + E.39 + E.40, 
                          data = LDP.OC.SNTV.pdata.Study1, 
                          model = "within")
round(summary(Study1.LDP.OC.SNTV, 
              vcov = vcovHC(Study1.LDP.OC.SNTV, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)
nrow(LDP.OC.SNTV.pdata.Study1)  # number of observations
length(unique(LDP.OC.SNTV.pdata.Study1$name_jp))  # number of MPs
length(unique(LDP.OC.SNTV.pdata.Study1$time))  # number of periods

## MMM, Budget Committee
Study1.LDP.BC.MMM <- plm(BC.char ~ leader + chair + 
                           minister + jr.minister + age + 
                           log.term + vote.share * SMD + zombie + 
                           E.42 + E.43 + E.44, 
                         data = LDP.BC.MMM.pdata.Study1, 
                         model = "within")
round(summary(Study1.LDP.BC.MMM, 
              vcov = vcovHC(Study1.LDP.BC.MMM, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)
nrow(LDP.BC.MMM.pdata.Study1)  # number of observations
length(unique(LDP.BC.MMM.pdata.Study1$name_jp))  # number of MPs
length(unique(LDP.BC.MMM.pdata.Study1$time))  # number of periods

round(exp(Study1.LDP.BC.MMM$coefficients[1]), 3)  # substantive interpretation

## MMM, other committees
Study1.LDP.OC.MMM <- plm(OC.char ~ leader + chair + 
                           minister + jr.minister + age + 
                           log.term + vote.share * SMD + zombie + 
                           E.42 + E.43 + E.44, 
                         data = LDP.OC.MMM.pdata.Study1, 
                         model = "within")
round(summary(Study1.LDP.OC.MMM, 
              vcov = vcovHC(Study1.LDP.OC.MMM, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)
nrow(LDP.OC.MMM.pdata.Study1)  # number of observations
length(unique(LDP.OC.MMM.pdata.Study1$name_jp))  # number of MPs
length(unique(LDP.OC.MMM.pdata.Study1$time))  # number of periods

#### robustness check: other dependent variables (Online Appendix D.1) ####
## number of speeches
Study1.LDP.PS.SNTV.count <- plm(PS.count ~ leader + chair + minister + 
                                  age + log.term + vote.share + 
                                  E.37 + E.38 + E.39 + E.40, 
                                data = LDP.PS.SNTV.pdata.Study1, 
                                model = "within")
round(summary(Study1.LDP.PS.SNTV.count, 
              vcov = vcovHC(Study1.LDP.PS.SNTV.count, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

Study1.LDP.BC.SNTV.count <- plm(BC.count ~ leader + chair + minister + 
                                  age + log.term + vote.share + 
                                  E.37 + E.38 + E.39 + E.40, 
                                data = LDP.BC.SNTV.pdata.Study1, 
                                model = "within")
round(summary(Study1.LDP.BC.SNTV.count, 
              vcov = vcovHC(Study1.LDP.BC.SNTV.count, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

Study1.LDP.OC.SNTV.count <- plm(OC.count ~ leader + chair + minister + 
                                  age + log.term + vote.share + 
                                  E.37 + E.38 + E.39 + E.40, 
                                data = LDP.OC.SNTV.pdata.Study1, 
                                model = "within")
round(summary(Study1.LDP.OC.SNTV.count, 
              vcov = vcovHC(Study1.LDP.OC.SNTV.count, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

Study1.LDP.PS.MMM.count <- plm(PS.count ~ leader + chair + 
                                 minister + jr.minister + age + 
                                 log.term + vote.share * SMD + zombie + 
                                 E.42 + E.43 + E.44, 
                               data = LDP.PS.MMM.pdata.Study1, 
                               model = "within")
round(summary(Study1.LDP.PS.MMM.count, 
              vcov = vcovHC(Study1.LDP.PS.MMM.count, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

Study1.LDP.BC.MMM.count <- plm(BC.count ~ leader + chair + 
                                 minister + jr.minister + age + 
                                 log.term + vote.share * SMD + zombie + 
                                 E.42 + E.43 + E.44, 
                               data = LDP.BC.MMM.pdata.Study1, 
                               model = "within")
round(summary(Study1.LDP.BC.MMM.count, 
              vcov = vcovHC(Study1.LDP.BC.MMM.count, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

Study1.LDP.OC.MMM.count <- plm(OC.count ~ leader + chair + 
                                 minister + jr.minister + age + 
                                 log.term + vote.share * SMD + zombie + 
                                 E.42 + E.43 + E.44, 
                               data = LDP.OC.MMM.pdata.Study1, 
                               model = "within")
round(summary(Study1.LDP.OC.MMM.count, 
              vcov = vcovHC(Study1.LDP.OC.MMM.count, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

## making one or more speeches
Study1.LDP.PS.SNTV.bin <- plm(PS.bin ~ leader + chair + minister + 
                                age + log.term + vote.share + 
                                E.37 + E.38 + E.39 + E.40, 
                              data = LDP.PS.SNTV.pdata.Study1, 
                              model = "within")
round(summary(Study1.LDP.PS.SNTV.bin, 
              vcov = vcovHC(Study1.LDP.PS.SNTV.bin, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

Study1.LDP.BC.SNTV.bin <- plm(BC.bin ~ leader + chair + minister + 
                                age + log.term + vote.share + 
                                E.37 + E.38 + E.39 + E.40, 
                              data = LDP.BC.SNTV.pdata.Study1, 
                              model = "within")
round(summary(Study1.LDP.BC.SNTV.bin, 
              vcov = vcovHC(Study1.LDP.BC.SNTV.bin, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

Study1.LDP.OC.SNTV.bin <- plm(OC.bin ~ leader + chair + minister + 
                                age + log.term + vote.share + 
                                E.37 + E.38 + E.39 + E.40, 
                              data = LDP.OC.SNTV.pdata.Study1, 
                              model = "within")
round(summary(Study1.LDP.OC.SNTV.bin, 
              vcov = vcovHC(Study1.LDP.OC.SNTV.bin, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

Study1.LDP.PS.MMM.bin <- plm(PS.bin ~ leader + chair + 
                               minister + jr.minister + age + 
                               log.term + vote.share * SMD + zombie + 
                               E.42 + E.43 + E.44, 
                             data = LDP.PS.MMM.pdata.Study1, 
                             model = "within")
round(summary(Study1.LDP.PS.MMM.bin, 
              vcov = vcovHC(Study1.LDP.PS.MMM.bin, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

Study1.LDP.BC.MMM.bin <- plm(BC.bin ~ leader + chair + 
                               minister + jr.minister + age + 
                               log.term + vote.share * SMD + zombie + 
                               E.42 + E.43 + E.44, 
                             data = LDP.BC.MMM.pdata.Study1, 
                             model = "within")
round(summary(Study1.LDP.BC.MMM.bin, 
              vcov = vcovHC(Study1.LDP.BC.MMM.bin, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

Study1.LDP.OC.MMM.bin <- plm(OC.bin ~ leader + chair + 
                               minister + jr.minister + age + 
                               log.term + vote.share * SMD + zombie + 
                               E.42 + E.43 + E.44, 
                             data = LDP.OC.MMM.pdata.Study1, 
                             model = "within")
round(summary(Study1.LDP.OC.MMM.bin, 
              vcov = vcovHC(Study1.LDP.OC.MMM.bin, 
                            method = "arellano", 
                            type = "HC2"))$coefficients, 3)

#### robustness check: censored regression (Online Appendix D.2) ####
# function to draw initial values
initial.fun.Tobit <- function(seed1, seed2, n.var, n.id) {
  set.seed(seed1)
  alpha <- rnorm(1)
  beta <- rnorm(n.var)
  gamma <- rnorm(n.id)
  sigma <- runif(1, 1, 2)
  tau <- runif(1, 1, 2)
  list(alpha = alpha, beta = beta, 
       gamma = gamma, sigma = sigma, tau = tau, 
       .RNG.name = "base::Wichmann-Hill", .RNG.seed = seed2)
}
# draw initial values
Study1.LDP.PS.SNTV.Tobit.inits <- list()
for (i in 1:3) {
  Study1.LDP.PS.SNTV.Tobit.inits[[i]] <- 
    initial.fun.Tobit(100 * i, 1000 * i, 12, 
                      max(as.numeric(factor(LDP.PS.SNTV.data.Study1$name_jp))))
}
# MCMC sampling
Study1.LDP.PS.SNTV.Tobit <- run.jags(model = "random_effects_Tobit_JAGS.R", 
                                     monitor = c("beta", "alpha", "sigma", "tau"), 
                                     data = list(y = ifelse(LDP.PS.SNTV.data.Study1$PS.char > 0, 
                                                            LDP.PS.SNTV.data.Study1$PS.char, NA), 
                                                 z = (LDP.PS.SNTV.data.Study1$PS.char > 0) * 1, 
                                                 x = scale(as.matrix(LDP.PS.SNTV.data.Study1[, c("leader", "chair", "minister", 
                                                                                                 "female", "age", "dynastic", 
                                                                                                 "log.term", "vote.share", 
                                                                                                 "E.37", "E.38", "E.39", "E.40")]), 
                                                           scale = FALSE), 
                                                 id = as.numeric(factor(LDP.PS.SNTV.data.Study1$name_jp)), 
                                                 N = nrow(LDP.PS.SNTV.data.Study1), 
                                                 n.var = 12, 
                                                 n.id = max(as.numeric(factor(LDP.PS.SNTV.data.Study1$name_jp)))), 
                                     inits = Study1.LDP.PS.SNTV.Tobit.inits, 
                                     n.chains = 3, burnin = 5000, sample = 5000, modules = "glm", 
                                     summarise = FALSE, thin = 1, method = "parallel")
# convergence check
print(which(gelman.diag(Study1.LDP.PS.SNTV.Tobit, multivariate = FALSE)[[1]][, 1] > 1.1))
# posterior summary
round(summary(Study1.LDP.PS.SNTV.Tobit), 3)
# posterior probability that the coefficient of party leader (beta[1]) is negative
round(mean(unlist(Study1.LDP.PS.SNTV.Tobit$mcmc[, 1]) < 0), 3)

Study1.LDP.BC.SNTV.Tobit.inits <- list()
for (i in 1:3) {
  Study1.LDP.BC.SNTV.Tobit.inits[[i]] <- 
    initial.fun.Tobit(100 * i, 1000 * i, 12, 
                      max(as.numeric(factor(LDP.BC.SNTV.data.Study1$name_jp))))
}
Study1.LDP.BC.SNTV.Tobit <- run.jags(model = "random_effects_Tobit_JAGS.R", 
                                     monitor = c("beta", "alpha", "sigma", "tau"), 
                                     data = list(y = ifelse(LDP.BC.SNTV.data.Study1$BC.char > 0, LDP.BC.SNTV.data.Study1$BC.char, NA), 
                                                 z = (LDP.BC.SNTV.data.Study1$BC.char > 0) * 1, 
                                                 x = scale(as.matrix(LDP.BC.SNTV.data.Study1[, c("leader", "chair", "minister", 
                                                                                                 "female", "age", "dynastic", 
                                                                                                 "log.term", "vote.share", 
                                                                                                 "E.37", "E.38", "E.39", "E.40")]), 
                                                           scale = FALSE), 
                                                 id = as.numeric(factor(LDP.BC.SNTV.data.Study1$name_jp)), 
                                                 N = nrow(LDP.BC.SNTV.data.Study1), 
                                                 n.var = 12, 
                                                 n.id = max(as.numeric(factor(LDP.BC.SNTV.data.Study1$name_jp)))), 
                                     inits = Study1.LDP.BC.SNTV.Tobit.inits, 
                                     n.chains = 3, burnin = 5000, sample = 5000, modules = "glm", 
                                     summarise = FALSE, thin = 1, method = "parallel")
print(which(gelman.diag(Study1.LDP.BC.SNTV.Tobit, multivariate = FALSE)[[1]][, 1] > 1.1))
round(summary(Study1.LDP.BC.SNTV.Tobit), 3)
round(mean(unlist(Study1.LDP.BC.SNTV.Tobit$mcmc[, 1]) < 0), 3)

Study1.LDP.OC.SNTV.Tobit.inits <- list()
for (i in 1:3) {
  Study1.LDP.OC.SNTV.Tobit.inits[[i]] <- 
    initial.fun.Tobit(100 * i, 1000 * i, 12, 
                      max(as.numeric(factor(LDP.OC.SNTV.data.Study1$name_jp))))
}
Study1.LDP.OC.SNTV.Tobit <- run.jags(model = "random_effects_Tobit_JAGS.R", 
                                     monitor = c("beta", "alpha", "sigma", "tau"), 
                                     data = list(y = ifelse(LDP.OC.SNTV.data.Study1$OC.char > 0, LDP.OC.SNTV.data.Study1$OC.char, NA), 
                                                 z = (LDP.OC.SNTV.data.Study1$OC.char > 0) * 1, 
                                                 x = scale(as.matrix(LDP.OC.SNTV.data.Study1[, c("leader", "chair", "minister", 
                                                                                                 "female", "age", "dynastic", 
                                                                                                 "log.term", "vote.share", 
                                                                                                 "E.37", "E.38", "E.39", "E.40")]), 
                                                           scale = FALSE), 
                                                 id = as.numeric(factor(LDP.OC.SNTV.data.Study1$name_jp)), 
                                                 N = nrow(LDP.OC.SNTV.data.Study1), 
                                                 n.var = 12, 
                                                 n.id = max(as.numeric(factor(LDP.OC.SNTV.data.Study1$name_jp)))), 
                                     inits = Study1.LDP.OC.SNTV.Tobit.inits, 
                                     n.chains = 3, burnin = 5000, sample = 5000, modules = "glm", 
                                     summarise = FALSE, thin = 1, method = "parallel")
print(which(gelman.diag(Study1.LDP.OC.SNTV.Tobit, multivariate = FALSE)[[1]][, 1] > 1.1))
round(summary(Study1.LDP.OC.SNTV.Tobit), 3)
round(mean(unlist(Study1.LDP.OC.SNTV.Tobit$mcmc[, 1]) < 0), 3)

Study1.LDP.PS.MMM.Tobit.inits <- list()
for (i in 1:3) {
  Study1.LDP.PS.MMM.Tobit.inits[[i]] <- 
    initial.fun.Tobit(100 * i, 1000 * i, 15, 
                      max(as.numeric(factor(LDP.PS.MMM.data.Study1$name_jp))))
}
Study1.LDP.PS.MMM.Tobit <- run.jags(model = "random_effects_Tobit_JAGS.R", 
                                    monitor = c("beta", "alpha", "sigma", "tau"), 
                                    data = list(y = ifelse(LDP.PS.MMM.data.Study1$PS.char > 0, LDP.PS.MMM.data.Study1$PS.char, NA), 
                                                z = (LDP.PS.MMM.data.Study1$PS.char > 0) * 1, 
                                                x = scale(as.matrix(LDP.PS.MMM.data.Study1[, c("leader", "chair", "minister", "jr.minister", 
                                                                                               "female", "age", "dynastic", 
                                                                                               "log.term", "vote.share", "SMD", "SMD.VS", "zombie", 
                                                                                               "E.42", "E.43", "E.44")]), 
                                                          scale = FALSE), 
                                                id = as.numeric(factor(LDP.PS.MMM.data.Study1$name_jp)), 
                                                N = nrow(LDP.PS.MMM.data.Study1), 
                                                n.var = 15, 
                                                n.id = max(as.numeric(factor(LDP.PS.MMM.data.Study1$name_jp)))), 
                                    inits = Study1.LDP.PS.MMM.Tobit.inits, 
                                    n.chains = 3, burnin = 5000, sample = 25000, modules = "glm", 
                                    summarise = FALSE, thin = 5, method = "parallel")
print(which(gelman.diag(Study1.LDP.PS.MMM.Tobit, multivariate = FALSE)[[1]][, 1] > 1.1))
round(summary(Study1.LDP.PS.MMM.Tobit), 3)
round(mean(unlist(Study1.LDP.PS.MMM.Tobit$mcmc[, 1]) < 0), 3)

Study1.LDP.BC.MMM.Tobit.inits <- list()
for (i in 1:3) {
  Study1.LDP.BC.MMM.Tobit.inits[[i]] <- 
    initial.fun.Tobit(100 * i, 1000 * i, 15, 
                      max(as.numeric(factor(LDP.BC.MMM.data.Study1$name_jp))))
}
Study1.LDP.BC.MMM.Tobit <- run.jags(model = "random_effects_Tobit_JAGS.R", 
                                    monitor = c("beta", "alpha", "sigma", "tau"), 
                                    data = list(y = ifelse(LDP.BC.MMM.data.Study1$BC.char > 0, LDP.BC.MMM.data.Study1$BC.char, NA), 
                                                z = (LDP.BC.MMM.data.Study1$BC.char > 0) * 1, 
                                                x = scale(as.matrix(LDP.BC.MMM.data.Study1[, c("leader", "chair", "minister", "jr.minister", 
                                                                                               "female", "age", "dynastic", 
                                                                                               "log.term", "vote.share", "SMD", "SMD.VS", "zombie", 
                                                                                               "E.42", "E.43", "E.44")]), 
                                                          scale = FALSE), 
                                                id = as.numeric(factor(LDP.BC.MMM.data.Study1$name_jp)), 
                                                N = nrow(LDP.BC.MMM.data.Study1), 
                                                n.var = 15, 
                                                n.id = max(as.numeric(factor(LDP.BC.MMM.data.Study1$name_jp)))), 
                                    inits = Study1.LDP.BC.MMM.Tobit.inits, 
                                    n.chains = 3, burnin = 5000, sample = 5000, modules = "glm", 
                                    summarise = FALSE, thin = 1, method = "parallel")
print(which(gelman.diag(Study1.LDP.BC.MMM.Tobit, multivariate = FALSE)[[1]][, 1] > 1.1))
round(summary(Study1.LDP.BC.MMM.Tobit), 3)
round(mean(unlist(Study1.LDP.BC.MMM.Tobit$mcmc[, 1]) < 0), 3)

Study1.LDP.OC.MMM.Tobit.inits <- list()
for (i in 1:3) {
  Study1.LDP.OC.MMM.Tobit.inits[[i]] <- 
    initial.fun.Tobit(100 * i, 1000 * i, 15, 
                      max(as.numeric(factor(LDP.OC.MMM.data.Study1$name_jp))))
}
Study1.LDP.OC.MMM.Tobit <- run.jags(model = "random_effects_Tobit_JAGS.R", 
                                    monitor = c("beta", "alpha", "sigma", "tau"), 
                                    data = list(y = ifelse(LDP.OC.MMM.data.Study1$OC.char > 0, LDP.OC.MMM.data.Study1$OC.char, NA), 
                                                z = (LDP.OC.MMM.data.Study1$OC.char > 0) * 1, 
                                                x = scale(as.matrix(LDP.OC.MMM.data.Study1[, c("leader", "chair", "minister", "jr.minister", 
                                                                                               "female", "age", "dynastic", 
                                                                                               "log.term", "vote.share", "SMD", "SMD.VS", "zombie", 
                                                                                               "E.42", "E.43", "E.44")]), 
                                                          scale = FALSE), 
                                                id = as.numeric(factor(LDP.OC.MMM.data.Study1$name_jp)), 
                                                N = nrow(LDP.OC.MMM.data.Study1), 
                                                n.var = 15, 
                                                n.id = max(as.numeric(factor(LDP.OC.MMM.data.Study1$name_jp)))), 
                                    inits = Study1.LDP.OC.MMM.Tobit.inits, 
                                    n.chains = 3, burnin = 5000, sample = 5000, modules = "glm", 
                                    summarise = FALSE, thin = 1, method = "parallel")
print(which(gelman.diag(Study1.LDP.OC.MMM.Tobit, multivariate = FALSE)[[1]][, 1] > 1.1))
round(summary(Study1.LDP.OC.MMM.Tobit), 3)
round(mean(unlist(Study1.LDP.OC.MMM.Tobit$mcmc[, 1]) < 0), 3)

#### robustness check: weighted linear fixed-effect model (Online Appendix D.2) ####
Study1.LDP.PS.SNTV.WFE <- wfe(PS.char ~ leader + chair + minister + 
                                age + log.term + vote.share + 
                                E.37 + E.38 + E.39 + E.40, 
                              data = LDP.PS.SNTV.data.Study1, 
                              unit.index = "name_jp", time.index = "time", 
                              treat = "leader", estimator = "fd")
round(summary(Study1.LDP.PS.SNTV.WFE)$coefficients, 3)

Study1.LDP.BC.SNTV.WFE <- wfe(BC.char ~ leader + chair + minister + 
                                age + log.term + vote.share + 
                                E.37 + E.38 + E.39 + E.40, 
                              data = LDP.BC.SNTV.data.Study1, 
                              unit.index = "name_jp", time.index = "time", 
                              treat = "leader", estimator = "fd")
round(summary(Study1.LDP.BC.SNTV.WFE)$coefficients, 3)

Study1.LDP.OC.SNTV.WFE <- wfe(OC.char ~ leader + chair + minister + 
                                age + log.term + vote.share + 
                                E.37 + E.38 + E.39 + E.40, 
                              data = LDP.OC.SNTV.data.Study1, 
                              unit.index = "name_jp", time.index = "time", 
                              treat = "leader", estimator = "fd")
round(summary(Study1.LDP.OC.SNTV.WFE)$coefficients, 3)

Study1.LDP.PS.MMM.WFE <- wfe(PS.char ~ leader + chair + 
                               minister + jr.minister + age + 
                               log.term + vote.share * SMD + zombie + 
                               E.42 + E.43 + E.44, 
                             data = LDP.PS.MMM.data.Study1, 
                             unit.index = "name_jp", time.index = "time", 
                             treat = "leader", estimator = "fd")
round(summary(Study1.LDP.PS.MMM.WFE)$coefficients, 3)

Study1.LDP.BC.MMM.WFE <- wfe(BC.char ~ leader + chair + 
                               minister + jr.minister + age + 
                               log.term + vote.share * SMD + zombie + 
                               E.42 + E.43 + E.44, 
                             data = LDP.BC.MMM.data.Study1, 
                             unit.index = "name_jp", time.index = "time", 
                             treat = "leader", estimator = "fd")
round(summary(Study1.LDP.BC.MMM.WFE)$coefficients, 3)

Study1.LDP.OC.MMM.WFE <- wfe(OC.char ~ leader + chair + 
                               minister + jr.minister + age + 
                               log.term + vote.share * SMD + zombie + 
                               E.42 + E.43 + E.44, 
                             data = LDP.OC.MMM.data.Study1, 
                             unit.index = "name_jp", time.index = "time", 
                             treat = "leader", estimator = "fd")
round(summary(Study1.LDP.OC.MMM.WFE)$coefficients, 3)